1.durée d’une simulation

1.1

ech<--log(runif(1000,0,1))
hist(ech,probability=T)
curve(dexp(x,1),col="darkgreen",add=T)

1.2

ech<-rexp(1000,1)
hist(ech,probability=T)
curve(dexp(x,1),col="darkgreen",add=T)

1.3

t0=proc.time()
a<-runif(10^6)
proc.time()-t0
##    user  system elapsed 
##   0.126   0.007   0.168

comparation

t0=proc.time()
ech<--log(runif(10000,0,1))
proc.time()-t0
##    user  system elapsed 
##   0.006   0.000   0.006
t0=proc.time()
ech<-rexp(10000,1)
proc.time()-t0
##    user  system elapsed 
##   0.005   0.000   0.006

2.Lien géométrique/exponentielle

2.1

Soit Yn suit une loi de géométrique, {Yn=k} ={B1=0,…,Bk=0,Bk+1=1}.

2.2

On veut simuler Bn qui suit une loi de Bernoulli.

p<-1/3
B1<-ceiling(runif(1)-(1-p))
count<-0
p<-1/3
while(B1==0){
  B1<-ceiling(runif(1)-(1-p))
  print(B1)
  count<-count+1
}
## [1] 1
print(count-1)
## [1] 0

count-1 suit une loi géométrique. Quand p est proche de 0, il prend bcp de fois avant réussir

2.3

Imx=N* ; P(x=k)=exp(-lambdak)-exp(-lambda(K+1))

on fait p=1-exp(-lambda)

t0=proc.time()
p=1/3
ech<-ceiling(log(runif(1000000))/log(1-p))
mean(ech)# théoriquement, il égale à 1/p
## [1] 3.003345
proc.time()-t0
##    user  system elapsed 
##   0.181   0.007   0.217
t0=proc.time()
p<-1/3
ech<-rgeom(1000000,p)
mean(ech)
## [1] 1.999768
proc.time()-t0
##    user  system elapsed 
##   0.219   0.004   0.246

On peut voir que ce qu’on fait est plus vite que rgeom de R! Pourquoi? parceque le R regard rgeom comme un cas spéciale de loi NB.

3.Nul n’est censé ignorer la loi normal

3.1

simulation d’une loi de Laplace

On va simule 100000 observations qui suit une loi normale. simuler une loi de bernoulli à valeur à {+1,-1}.

simuler une loi de Laplace

m<-sqrt(2*exp(1)/pi)
r<- function (x) exp(-0.5*((abs(x)-1)^2))
c<-c()
t0=proc.time()
count<-0
while(count<10000){ 
ech_ber<-(ceiling((runif(1)-1/2))-1/2)*2
ech<--log(runif(1))*ech_ber
if (runif(1)<r(ech)) {
  c<-append(c,ech)
  count<-count+1
}}
proc.time()-t0
##    user  system elapsed 
##   0.743   0.051   0.954
hist(c,probability=T)
curve(dnorm(x),add=T,col="red")
lines(density(c),col="darkgreen")
legend("topright",c("desité théprique","densité empirique"),col=c("red","darkgreen"),lty=c(1,1))

Ils coincident parfaitement.

3.2

la méthode de Box-Muller.

t0<-proc.time()
n<-50000
r<-sqrt(-2*log(runif(n)))
theta<-2*pi*runif(n)
N2<-c(r*cos(theta),r*sin(theta))
hist(N2,probability=T)
curve(dnorm(x),add=T,col="red")
lines(density(N2),col="darkgreen")
legend("topright",c("desité théprique","densité empirique"),col=c("red","darkgreen"),lty=c(1,1))

proc.time()-t0
##    user  system elapsed 
##   0.082   0.009   0.144

Cette méthode est bcp plus vite que la méthode rejet que l’on a fait.

3.3

n<-1000
r<-sqrt(-2*log(runif(n)))
theta<-2*pi*runif(n)
x1<-r*cos(theta)
x2<-r*sin(theta)
X<-rbind(x1,x2)

X est une échantillon qui suit une loi normale std de dim 2.

gamma=matrix(c(1,1,1,4),nrow=2,byrow = T)
library(expm)#le package est utilisé pour calculer le sqrt(resp. log, exp) d'une matice.
## Loading required package: Matrix
## 
## Attaching package: 'expm'
## 
## The following object is masked from 'package:Matrix':
## 
##     expm
A<-sqrtm(gamma)
## Note: method with signature 'symmetricMatrix#missing' chosen for function 'Schur',
##  target signature 'dsyMatrix#missing'.
##  "dsyMatrix#ANY" would also be valid
B<-matrix(c(2,1),nrow = 2)
f<-function (x) A%*%x+B
Y<-apply(X,2,f)
plot(Y[1,],Y[2,],col="darkred",xlim=c(-2,5),ylim=c(-6,8),main="Représentation de Y",
     xlab="Y1",ylab="Y2")

On voit bien une ellipsoid.

apply(Y,1,mean)
## [1] 1.9597466 0.9948785
cov(t(Y))
##          [,1]     [,2]
## [1,] 1.000922 1.109956
## [2,] 1.109956 4.092490

Si on veut utiliser directement le package “MASS”

n=1000
m=c(2,1)
gamma=matrix(c(1,1,1,4),nrow=2)
library(MASS)
ech=mvrnorm(n,m,gamma)
cov(ech)
##          [,1]     [,2]
## [1,] 1.013282 1.085869
## [2,] 1.085869 4.281866
plot(ech[,1],ech[,2],col="darkred",xlim=c(-2,5),ylim=c(-6,8),main="Représentation de ech",
     xlab="ech1",ylab="ech2")

4.Méfiez-vous des mélanges !

4.1

f<- function (x) (1/(6*sqrt(2*pi)))*exp(-(x+3)^2/8)+(2/(3*sqrt(2*pi))*exp(-(x-3)^2/2))
plot(f,xlim = c(-10,10),main="Représentation de f")

4.2&4.3

on regard X comme une loi mélange de N(-3,2) et N(3,1)

u<-runif(10000)
X<-(u<1/3)*rnorm(10000,-3,2)+(u>1/3)*rnorm(10000,3,1)
hist(X,probability = T)
lines(density(X),col="blue")
curve(f(x),add=T,col="red")
legend("topleft",c("densité empirique","densité vraie"),col = c("blue","red"),lty = c(1,1))

Il coincide parfaitement !

5.Représentation et simulation de copules.

5.1

Initiation

library(rgl)
x<-seq(0.01,1,0.01)
y<-x

copule_min

c_min<- function(x,y) {
  X<-rbind(x,y)
  return(apply(X,2,min))
}
c_min_z<-outer(x,y,c_min)
#surface3d(x, y, c_min_z,col="red")
persp(x,y,c_min_z,col="red")

copule anti-comonotone

c_max<- function(x,y) {
  X<-rbind(x,y)
  f<-function(u){ max(u[1]+u[2]-1,0)}
  return(apply(X,2,f))
}
c_max_z<-outer(x,y,c_max)
persp(x,y,c_max_z,col="red")

#surface3d(x, y, c_max_z,col="red")

copule indépendante

c_indep_z<-outer(x,y,"*")
persp(x,y,c_indep_z,col="red")

#surface3d(x, y, c_indep_z,col="red")

copule de clayton(de paramètre 1)

theta=1
c_clayton<- function(x,y) {
  X<-rbind(x,y)
  f<-function(u){ (u[1]^(-theta)+u[2]^(-theta)-1)^(-1/theta)}
  return(apply(X,2,f))
}

Attention! x et y ne peuvent pas prendre 0

c_clayton_z<-outer(x,y,c_clayton)
persp(x,y,c_clayton_z,col="red")

#surface3d(x, y, c_clayton_z,col="red")

très joli!

copule gaussien.

p=1/2
library(mvtnorm)
library(rgl)
x<-seq(0.01,1,0.01)
y<-x
c_gaussien<-function(x,y) {
  X<-rbind(x,y)
  f<-function(u){ 
    mean <- rep(0, 2)
    lower<- -Inf
    upper <- c(qnorm(min(1,u[1])),qnorm(min(1,u[2])))
    corr <- matrix(c(1,p,p,1),nrow = 2)
    pmvnorm(lower,upper,mean,corr) }
  return(apply(X,2,f))
}

c_gaussien_z<-outer(x,y,c_gaussien)
persp(x,y,c_gaussien_z,col="red")

#surface3d(x, y, c_gaussien_z,col="red")

5.2

library(MASS)
p<-1/2
gamma <- matrix(c(1,p,p,1),nrow = 2)
X=mvrnorm(1000,rep(0,2),gamma)
u<-pnorm(X[,1])
v<-pnorm(X[,2])

(u,v)est une échantiollon dont la fonction de repartion égale à la copule gaussien de paramètre p.

U<-runif(1000)
u<-runif(1000)
V<-U*sqrt(u)/(1+U*sqrt(u)-sqrt(u))

(U,V)est un couple de uniform dont la fdr est c_clayton.

par(mfrow=c(1,2))
plot(X,main = "Copule Gaussien",col='darkblue')
plot(qnorm(U),qnorm(V),main = 'Copule de Clayton',col='darkgreen')

On peut voir qu’ils se ressemblent un peu.

5.3

Simulation d’une loi de Laplace

n<-1000
ech_ber1<-2*(ceiling(runif(n)-1/2)-1/2)
ech_exp1<--log(runif(n))
ech_lap1<-ech_ber1*ech_exp1
ech_ber2<-2*(ceiling(runif(n)-1/2)-1/2)
ech_exp2<--log(runif(n))
ech_lap2<-ech_ber2*ech_exp2
par(mfrow=c(1,2))

copule indépendante & copule gaussienne

library(MASS)
p<-1/2
gamma <- matrix(c(1,p,p,1),nrow = 2)
X=mvrnorm(1000,rep(0,2),gamma)
u<-pnorm(X[,1])
v<-pnorm(X[,2])

(u,v)est une échantiollon dont la fonction de repartion égale à la copule gaussien de paramètre p.

x<-rep(0,n)
y<-rep(0,n)
for(i in 1:n){
  if (u[i]<1/2){
    x[i]<-log(2*u[i])
  }else{
    x[i]<--log(2*(1-u[i]))
}}
for(i in 1:n){
  if (v[i]<1/2){
    y[i]<-log(2*v[i])
  }else{
    y[i]<--log(2*(1-v[i]))
  }}
par(mfrow=c(1,2))
plot(ech_lap1,ech_lap2,main = 'Copule Indépendante',col='darkblue')
plot(x,y,main = 'Copule Gussienne',col='darkgreen')

6.Movement brownien, pont brownien.

6.1

n<-10000
t<-seq(1/n,1,1/n)
u<-rnorm(n,0,sqrt(1/n))
W<-cumsum(u)
plot(t,W,type = 'l',col="darkred",main = 'Movement de Brownien')

6.2

n<-10000
u1<-rnorm(n,0,sqrt(1/n))
W1<-cumsum(u1)
u2<-rnorm(n,0,sqrt(1/n))
W2<-cumsum(u2)
plot(W1,W2,type = 'l',col="darkgreen",main = 'Movement de Brownien(dim 2)')
points(0,0,col="red")

6.3

n<-10000
t<-seq(1/n,1,1/n)
u<-rnorm(n,0,sqrt(1/n))
W<-cumsum(u)
B<-rep(0,n)
for( i in 1:n) B[i]<-W[i]-t[i]*W[n]
plot(t,B,type = 'l',col="darkred",main = 'Pont de Brownien')

6.4(a)

n=10000
U<-runif(n)
U<-sort(U)
Fn<- function(x) 1/n*sum((rep(x,n)-U>=0))
x=seq(1/n,1,by=1/n)
x<-matrix(x,nrow=1)
plot(x,apply(x,2,Fn),type ='s')

FB<-function(x) sqrt(n)*(Fn(x)-x)#c'est la fonction sqrt(n)*(Fn(x)-F(x))
plot(x,apply(x,2,FB),type ='s')

6.4(b)

n<-100
x=seq(1/n,1,by=1/n)
x<-matrix(x,nrow=1)
#On simule une échantillons de loi sqrt(n)*||Fn-F||inf de taille 1000
X_B<-c()
FB_abs<-function(x) sqrt(n)*abs((Fn(x)-x))
for (i in 1:1000){
  U<-runif(n)
  U<-sort(U)
  Fn<- function(x) 1/n*sum((rep(x,n)-U>=0))
  X_B<-append(X_B,max(apply(x,2,FB)))
}

On simule une échantillons de loi sup(Bt) de taille 1000

sup_B=c()
for(i in 1:1000){
  t<-seq(1/n,1,1/n)
  u<-rnorm(n,0,sqrt(1/n))
  W<-cumsum(u)
  B<-rep(0,n)
  for( i in 1:n) B[i]<-W[i]-t[i]*W[n]
  sup_B<-append(sup_B,max(B))
}

comparaison!

par(mfrow=c(1,2))
hist(X_B,probability = T)
lines(density(X_B),col='red')
lines(density(sup_B),col='blue')
legend("topright",c("X_B","sup_B"),col=c("red","blue"),lty=c(1,1))

hist(sup_B,probability = T)
lines(density(X_B),col='red')
lines(density(sup_B),col='blue')
legend("topright",c("X_B","sup_B"),col=c("red","blue"),lty=c(1,1))

C’est magic !